Predicting Aquaculture Suitability

Author

Sam Muir

Published

December 4, 2023

Overview

Marine aquaculture has the potential to play an important role in the global food supply as a more sustainable protein option than land-based meat production.1 Gentry et al. mapped the potential for marine aquaculture globally based on multiple constraints, including ship traffic, dissolved oxygen, bottom depth .2

In this project, I am working to determine which Exclusive Economic Zones (EEZ) on the West Coast of the US are best suited to developing marine aquaculture for several species of oysters.

Based on previous research, we know that oysters needs the following conditions for optimal growth:

  • sea surface temperature: 11-30°C
  • depth: 0-70 meters below sea level

From this analysis, I will also develop a function to determine which Exclusive Economic Zones (EEZ) on the West Coast of the US are best suited for other aquaculture species based on their temperature and depth ranges.

Steps of Analysis:
  • combining vector/raster data
  • resampling raster data
  • masking raster data
  • map data
  • generalize to a reusable function

Data

Sea Surface Temperature

We will use average annual sea surface temperature (SST) from the years 2008 to 2012 to characterize the average sea surface temperature within the region. The data we are working with was originally generated from NOAA’s 5km Daily Global Satellite Sea Surface Temperature Anomaly v3.1.

Bathymetry

To characterize the depth of the ocean we will use the General Bathymetric Chart of the Oceans (GEBCO).3

Exclusive Economic Zones

We will be designating maritime boundaries using Exclusive Economic Zones off of the west coast of US from Marineregions.org.

Prepare data

To start, we need to load all necessary data and make sure they have matching coordinate reference systems.

# load libraries
library(tidyverse)
library(here)
library(sf)
library(terra)
library(ggspatial)
library(leaflet)
library(htmltools)
# read in shapefiles
eez <- read_sf(here("data", "wc_regions_clean.shp"))

# read in SST rasters
sst_2008 <- rast(here("data", "average_annual_sst_2008.tif"))
sst_2009 <- rast(here("data", "average_annual_sst_2009.tif"))
sst_2010 <- rast(here("data", "average_annual_sst_2010.tif"))
sst_2011 <- rast(here("data", "average_annual_sst_2011.tif"))
sst_2012 <- rast(here("data", "average_annual_sst_2012.tif"))

# combining the sst data 
sst <- list(sst_2008, sst_2009, sst_2010, sst_2011, sst_2012)
sst_rast <- rast(sst)
# plotting check
#plot(sst_rast)

# read in bathymetry raster
depth <- rast(here("data", "depth.tif"))
# plotting check
#plot(depth)

# getting base map for the US/west coast
state <- map_data('state')
westcoast <- subset(state, region %in% c("washington", "oregon", "california"))
# check crs of all data
#st_crs(wc) # 4326
#st_crs(sst_rast) # 9122
#st_crs(depth) # 4326

# reproject sst_rast to have same crs as the others
sst_rast <- project(sst_rast, "EPSG:4326")

# check to see that the crs matches
# st_crs(sst_rast) == st_crs(depth)

Process data

Next, we need process the SST and depth data so that they can be combined. In this case the SST and depth data have slightly different resolutions, extents, and positions. We don’t want to change the underlying depth data, so we will need to resample to match the SST data using the nearest neighbor approach.

# mean sst from 2008-2012
sst_mean <-  mean(sst_rast)

# calculating the mean K to ºC
sst_mean_c <- sst_mean - 273.15

# cropping depth data to sst raster extent
depth_crop <- crop(depth, sst_mean_c)

# matching the resolution of both rasters
depth_crop_res <- resample(x = depth_crop, 
                               y = sst_mean_c, 
                               method = "near")

#checking that the CRS, resolution, and extent match
st_crs(sst_mean_c) == st_crs(depth_crop_res)
[1] TRUE
res(sst_mean_c) == res(depth_crop_res)
[1] TRUE TRUE
ext(sst_mean_c) == ext(depth_crop_res)
[1] TRUE
# stacking the two rasters of depth and sea surface temperature
sst_depth_list <- list(sst_mean_c, depth_crop_res)
sst_depth_rast <- rast(sst_depth_list)

# plotting new raster
# plot(sst_depth_rast)

Find suitable locations

In order to find suitable locations for marine aquaculture, we’ll need to find locations that are suitable in terms of both SST and depth.

# creating a reclassification matrix valid sst locations
sst_reclass <- matrix(c(-Inf, 11, NA, 
                        11, 30, 1,
                        30, Inf, NA),
                      ncol = 3,
                      byrow = TRUE)

# using the reclassifying matrix to set non-suitable sst to NA
sst_suitable <- classify(sst_mean_c, rcl = sst_reclass)

# creating a reclassification matrix valid depth locations
depth_reclass <- matrix(c(-Inf, -70, NA, 
                          -70, 0, 1,
                          0, Inf, NA),
                        ncol = 3,
                        byrow = TRUE)

# using the reclassifying matrix to set non-suitable depth to NA
depth_suitable <- classify(depth_crop_res, rcl = depth_reclass)

# cropping sst and depth data based on the mask
sst_oyster <- mask(sst_mean_c, sst_suitable)
depth_oyster <- mask(depth_crop_res, depth_suitable)

# combining the two rasters
list_oyster <- list(depth_suitable, sst_suitable)
sst_depth_oyster <- rast(list_oyster)

# overlaying the data
fun_oyster <- function(x, y){
  return(x * y)
  }

sst_depth_suitable <- lapp(sst_depth_oyster,
                               fun = fun_oyster)

# check by plotting
#plot(sst_depth_suitable)

Determine the most suitable EEZ

We want to determine the total suitable area within each EEZ in order to rank zones by priority. To do so, we need to find the total area of suitable locations within each EEZ.

#mask cell size for suitable areas
mask_suitable <- cellSize(sst_depth_suitable,
                          mask = TRUE,
                          unit = 'km',
                          transform = T)

#rasterize wc region
wc_rasterized <- rasterize(eez, sst_depth_suitable, field = 'rgn')

#total suitable area
total_suitable_area <- zonal(mask_suitable, wc_rasterized,
                            fun = 'sum',
                            na.rm = TRUE)

#percent suitable area
percentage_suitable_area <- left_join(eez, total_suitable_area, by = 'rgn') %>%
  mutate(percent_suitable_area = (area / area_km2) * 100)

# report in a table
percentage_suitable_area %>%
  st_drop_geometry() %>%
  select(Region = rgn,
         `Total Suitable Area` = area_km2,
         `Percent Suitable Area` = percent_suitable_area)
# A tibble: 5 × 3
  Region              `Total Suitable Area` `Percent Suitable Area`
  <chr>                               <dbl>                   <dbl>
1 Oregon                            179994.                   0.597
2 Northern California               164379.                   0.108
3 Central California                202738.                   2.01 
4 Southern California               206861.                   1.70 
5 Washington                         66898.                   3.56

Plotting the results

pal <- colorBin("Blues", domain = percentage_suitable_area$percent_suitable_area) # define color palette
pal_total <- colorBin("Blues", domain = percentage_suitable_area$area_km2)
# labels for the hover over text
labels_percent <- sprintf(
  "<strong>%s</strong><br/>%g percent",
  percentage_suitable_area$rgn, round(percentage_suitable_area$percent_suitable_area, 2)) %>% 
  lapply(htmltools::HTML)
labels_total <- sprintf(
  "<strong>%s</strong><br/>%g km<sup>2</sup>",
  percentage_suitable_area$rgn, percentage_suitable_area$area_km2) %>% 
  lapply(htmltools::HTML)

# percent area map
perc <- leaflet(percentage_suitable_area) %>%
  addPolygons(
    fillColor = ~pal(percentage_suitable_area$percent_suitable_area), # fill by the percent suitable area
    weight = 2,
    opacity = 1,
    color = "grey", # border
    dashArray = "3",
    fillOpacity = 0.7,
    highlightOptions = highlightOptions( # highlight region when hovered over
    weight = 5,
    color = "#666",
    dashArray = "",
    fillOpacity = 0.7,
    bringToFront = TRUE),
  label = labels_percent,
  labelOptions = labelOptions( # labels when hover
    style = list("font-weight" = "normal", padding = "3px 8px"),
    textsize = "15px",
    direction = "auto")) %>%
  addLegend(pal = pal, # legend based on color palette
            values = ~percent_suitable_area, 
            opacity = 0.7,
            labFormat = labelFormat(suffix = "%"), # add units
            title = "% Suitable Area <br> for Oysters",
            position = "bottomright") %>%
  addTiles() # basemap

# total area map
total <- leaflet(percentage_suitable_area) %>%
  addPolygons(
    fillColor = ~pal_total(percentage_suitable_area$area_km2), # fill by the total suitable area
    weight = 2,
    opacity = 1,
    color = "grey",
    dashArray = "3",
    fillOpacity = 0.7,
   # popup = p_popup,
    highlightOptions = highlightOptions( # highlight when hoverover
    weight = 5,
    color = "#666",
    dashArray = "",
    fillOpacity = 0.7,
    bringToFront = TRUE),
  label = labels_total, # hover labels
  labelOptions = labelOptions(
    style = list("font-weight" = "normal", padding = "3px 8px"),
    textsize = "15px",
    direction = "auto")) %>%
  addLegend(pal = pal_total, # add legend with color palette
            values = ~area_km2, 
            opacity = 0.7,
            labFormat = labelFormat(suffix = " km<sup>2</sup>"), # add units
            title = "Total Suitable Area <br> for Oysters",
            position = "bottomright") %>%
  addTiles() # basemap

# put maps side by side
oyster_grid <- 
  tagList(
    tags$table(width = "100%",
      tags$tr(
        tags$td(perc),
        tags$td(total)
      )
      )
    )
# show maps
browsable(oyster_grid)

Creating a Re-Usable Function

species_suitability <- function(min_temp, max_temp, min_depth, max_depth, species) {

# Stack the rasters
sst <- list(sst_2008, sst_2009, sst_2010, sst_2011, sst_2012)
sst_rast <- rast(sst)

# transform crs
eez <- st_transform(eez, "EPSG:4326")
sst_rast <- project(sst_rast, "EPSG:4326")
depth <- project(depth, "EPSG:4326")

# find mean and convert to ºC
mean_sst <- mean(sst_rast)
mean_sst <- mean_sst - 273.15

#crop and resample
depth_crop <- crop(depth, mean_sst)
depth_resamp <- resample(x = depth_crop, y = mean_sst, method = "near")
#Reclassify depth raster to = 1 when between -70 and 0
depth_reclass <- matrix(c(-Inf, -max_depth, NA,
                      -max_depth, -min_depth, 1,
                      -min_depth, Inf, NA), ncol = 3, byrow = TRUE)
suitable_depth <- classify(depth_resamp, rcl = depth_reclass)
# Reclassify sst raster using funtion input of temps and depths
sst_reclass <- matrix(c(-Inf, min_temp, NA,
                    min_temp, max_temp, 1,
                    max_temp, Inf, NA), ncol = 3, byrow = TRUE)
suitable_sst <- classify(mean_sst, rcl = sst_reclass)
# Find suitable locations for both depth and sst
suitable_stack <- c(suitable_depth, suitable_sst)
fun_suitable <- function(x, y) {x*y}
suitable <- lapp(suitable_stack[[c(1, 2)]], fun = fun_suitable)
# mask cell size for suitable areas
mask_suitable <- cellSize(suitable,
                          mask = TRUE,
                          unit = 'km',
                          transform = T)

#rasterize eez
wc_rasterized <- rasterize(eez, suitable, field = 'rgn')

#total suitable area
total_suitable_area <- zonal(mask_suitable, wc_rasterized,
                            fun = 'sum',
                            na.rm = TRUE)

#percent suitable area
suitable_area <- left_join(eez, total_suitable_area, by = 'rgn') %>%
  mutate(percent_suitable_area = (area / area_km2) * 100)

# table
print(suitable_area %>%
  st_drop_geometry() %>%
  select(Region = rgn,
         `Percent Suitable Area` = percent_suitable_area,
         `Total Suitable Area` = area_km2))

# mapping
pal <- colorBin("Blues", domain = suitable_area$percent_suitable_area) # define color palette
pal_total <- colorBin("Blues", domain = suitable_area$area_km2)
# make labels for hover over
labels_percent <- sprintf(
  "<strong>%s</strong><br/>%g percent",
  suitable_area$rgn, round(suitable_area$percent_suitable_area, 2)) %>% 
  lapply(htmltools::HTML)
labels_total <- sprintf(
  "<strong>%s</strong><br/>%g km<sup>2</sup>",
  suitable_area$rgn, suitable_area$area_km2) %>% 
  lapply(htmltools::HTML)


percent <- leaflet(suitable_area) %>%
  addPolygons(
    fillColor = ~pal(suitable_area$percent_suitable_area), # color region by percent
    weight = 2,
    opacity = 1,
    color = "grey",
    dashArray = "3",
    fillOpacity = 0.7,
   # popup = p_popup,
    highlightOptions = highlightOptions( # when you hover over there will be pop up text
    weight = 5,
    color = "#666",
    dashArray = "",
    fillOpacity = 0.7,
    bringToFront = TRUE),
  label = labels_percent, # label by previously defined
  labelOptions = labelOptions(
    style = list("font-weight" = "normal", padding = "3px 8px"),
    textsize = "15px",
    direction = "auto")) %>%
  addLegend(pal = pal, # legend based on color palette
            values = ~suitable_area, 
            opacity = 0.7,
            labFormat = labelFormat(suffix = "%"), # add units
            title = paste("% Suitable Area <br> for", species), # updates with species input
            position = "bottomright") %>%
  addTiles() # basemap

total_area <- leaflet(suitable_area) %>%
  addPolygons(
    fillColor = ~pal_total(suitable_area$area_km2), # color region by area
    weight = 2,
    opacity = 1,
    color = "grey",
    dashArray = "3",
    fillOpacity = 0.7,
    highlightOptions = highlightOptions( # when you hover over there will be pop up text
    weight = 5,
    color = "#666",
    dashArray = "",
    fillOpacity = 0.7,
    bringToFront = TRUE),
  label = labels_total, # label by previously defined
  labelOptions = labelOptions(
    style = list("font-weight" = "normal", padding = "3px 8px"), # styling
    textsize = "15px",
    direction = "auto")) %>%
  addLegend(pal = pal_total, # legend based on color palette
            values = ~area_km2, 
            opacity = 0.7,
            labFormat = labelFormat(suffix = " km<sup>2</sup>"), # add units
            title = paste("Total Suitable Area <br> for", species), # title that updates with species input
            position = "bottomright") %>%
  addTiles() # basemap

# put the leaflet maps side by side
aquaculture_grid <- 
  tagList(
    tags$table(width = "100%",
      tags$tr(
        tags$td(percent),
        tags$td(total_area)
      )
      )
    )
# show the maps
browsable(aquaculture_grid)

}

Run the function with your aquaculture species of choice

SeaLifeBase is a good resource for finding information on species depth and temperature requirements.

# give the function a try!
species_suitability(min_temp = 10, max_temp = 14, 
                   min_depth =  50, max_depth = 500, 
                   species = "Blue Squat Lobster")
# A tibble: 5 × 3
  Region              `Percent Suitable Area` `Total Suitable Area`
  <chr>                                 <dbl>                 <dbl>
1 Oregon                                11.7                179994.
2 Northern California                    5.20               164379.
3 Central California                     5.14               202738.
4 Southern California                    1.54               206861.
5 Washington                            16.4                 66898.

Footnotes

  1. Hall, S. J., Delaporte, A., Phillips, M. J., Beveridge, M. & O’Keefe, M. Blue Frontiers: Managing the Environmental Costs of Aquaculture (The WorldFish Center, Penang, Malaysia, 2011).↩︎

  2. Gentry, R. R., Froehlich, H. E., Grimm, D., Kareiva, P., Parke, M., Rust, M., Gaines, S. D., & Halpern, B. S. Mapping the global potential for marine aquaculture. Nature Ecology & Evolution, 1, 1317-1324 (2017).↩︎

  3. GEBCO Compilation Group (2022) GEBCO_2022 Grid (doi:10.5285/e0f0bb80-ab44-2739-e053-6c86abc0289c).↩︎